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EFFICIENT PDE CONSTRAINED SHAPE OPTIMIZATION BASED 
ON STEKLOV-POINCARE TYPE METRICS 


VOLKER H. SCHULZ, MARTIN SIEBENBORN AND KATHRIN WELKER* 

Abstract. Recent progress in PDE constrained optimization on shape manifolds is based 
on the Hadamard form of shape derivatives, i.e., in the form of integrals at the boundary of the 
shape under investigation, as well as on intrinsic shape metrics. From a numerical point of view, 
domain integral forms of shape derivatives seem promising, which rather require an outer metric 
on the domain surrounding the shape boundary. This paper tries to harmonize both points of 
view by employing a Steklov-Poincare type intrinsic metric, which is derived from an outer metric. 
Based on this metric, efficient shape optimization algorithms are proposed, which also reduce the 
analytical labor, so far involved in the derivation of shape derivatives. 
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1. Introduction. Shape optimization is of interest in many fields of appli¬ 
cation - in particular in the context of partial differential equations (PDE). As 
examples, we mention aerodynamic shape optimization [22]’ acoustic shape opti¬ 
mization 130] or optimization of interfaces in transmission problems mmm and 
in electrostatics [2]. In industry, shapes are often represented within a finite dimen¬ 
sional design space. However, often this reduction is felt as being too restrictive 
mi which motivates shape optimization based on shape calculus. Major effort in 
shape calculus 0126] has been devoted towards expressions for shape derivatives 
in so-called Hadamard-form, i.e., in boundary integral form. It is known that the 
second order shape derivative, formerly coined as shape Hessian, is nonsymmetric in 
general, which for a long time has been an obstacle for algorithmic developments in 
shape optimization in the fashion of nonlinear programming. Recently |231 1241 123] , 
shape optimization has been considered as optimization on Riemannian shape man¬ 
ifolds, which enables design and analysis of NLP-like algorithms including one-shot 
sequential quadratic programming and theoretical insights into the structure of the 
second order shape derivative in comparison to the Riemannian shape Hessian. Co- 
ercivity results for shape Hessians for elliptic problems can be found in 0- The 
scalar product used in this work is in line with these results. 

On the other hand, it is often a very tedious, not to say painful, process to 
derive the boundary formulation of the shape derivative. Along the way, there fre¬ 
quently appears a domain formulation in the form of an integral over the whole 
domain as an intermediate step. Recently, it has been shown that this intermediate 
formulation has numerical advantages 0 noma [20]. In [13], also practical advan¬ 
tages of the domain shape formulation have been demonstrated, since it requires 
less smoothness assumptions. Furthermore, the derivation as well as the implemen¬ 
tation of the domain integral formulation requires less manual and programming 
work. Thus, there arises the natural goal of combining the favorable domain in¬ 
tegral formulation of shape derivatives with the favorable NLP-type optimization 
strategies on shape manifolds, which seem so far tightly coupled with boundary 
integral formulations of shape derivatives. This publication aims at demonstrating 
that this coupling is indeed possible and that it naturally leads to a novel family 
of Poincare-Steklov type metrics on shape manifolds. In contrast to [24] this work 
consciously avoids surface formulations of shape derivatives in order to provide more 
handy optimization algorithms. 

The paper is organized in the following way. First, in section [2] we set up 
notation and terminology and formulate the model problem. In section [3j we dis¬ 
cuss generalized Poincare-Stcklov operators as the basis for Riemannian metrics for 
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shape manifolds. Section [6] is devoted to the set of all shapes in the context of the 
novel metric introduced in section [3] Section [4] rephrases NLP-like optimization al¬ 
gorithms on shape manifolds within the framework of domain integral formulations 
of shape derivatives. Finally, section [5] discusses algorithmic and implementation 
details, as well as, numerical results for a parabolic transmission shape optimization 
problems. 

2. Problem Formulation. We first set up notation and terminology in shape 
calculus. Then we recall the model problem in |24] . which is motivated by electri¬ 
cal impedance tomography and given by a parabolic interface shape optimization 
problem. 

2.1. Notations in shape calculus. Let d £ N and r > 0. We denote by 
fl C a bounded domain with Lipschitz boundary T := dfl and by J a real¬ 
valued functional depending on it. Moreover, let {-F)}te[o,T] be a family of bijective 
mappings F t : fl —> such that Fq = id. This family transforms the domain fl 
into new perturbed domains fit '■= F t (fl) = { F t (x ): x £ f 1} with flo = fl and the 
boundary T into new perturbed boundaries T t := F t (T) = {F t (x): x £ T} with 
T 0 = T. If you consider the domain Q as a collection of material particles, which 
are changing their position in the time-interval [0,r], then the family {F)} tg [ 0T ] 
describes the motion of each particle, i.e., at the time t £ [0,r] a material particle 
x £ fl has the new position xt := F t (x) £ f l t with xq = x. The motion of each such 
particle x could be described by the velocity method , i.e., as the flow F t (x) := £(f, x) 
determined by the initial value problem 


£(0, x) = X 


( 2 . 1 ) 


or by the perturbation of identity , which is defined by F t (x) := x + tV(x) where V 
denotes a sufficiently smooth vector field. We will use the perturbation of identity 
throughout the paper. The Eulerian derivative of J at Q in direction V is defined 
by 


DJ(fl)\V] := lim 

t-> o+ 


j(n t ) - j{n) 

t 


( 2 . 2 ) 


The expression DJ(fl)[V] is called the shape derivative of J at fl in direction V 
and J shape differentiable at f l if for all directions V the Eulerian derivative (2.2) 
exists and the mapping V DJ(fl)[V ] is linear and continuous. For a thorough 
introduction into shape calculus, we refer to the monographs [7] 126] . In particular, 
m states that shape derivatives can always be expressed as boundary integrals due 
to the Hadamard structure theorem. The shape derivative arises in two equivalent 
notational forms: 

DJn[V]:= j F{x)V{x)dx (domain formulation) (2.3) 

Jfi 

DJr[V] := J f (s)V (s) T n(s) ds (boundary formulation) (2.4) 


where F(x) is a (differential) operator acting linearly on the perturbation vector 
field V and /: F —► R. with 


DJ n [V] = DJ(fL) [V] = DJ r [V). (2.5) 

The boundary formulation ( |2.4| ), DJr[V], acting on the normal component of V 
has led to the interpretation as tangential vector of a corresponding shape manifold 
in |23]. 
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2.2. PDE model definition. We use the same model problem as in Til , 
which is briefly recalled. Let this domain 12 be an open subset of R 2 and split into 
the two disjoint subdomains 12 i, 122 C 12 such that 12f U U 122 = 12, Tbottom LI 
rieftUr ri g ht ur t op = <912 (=: r out ) and <912in9122 = rj nt where the interior boundary 
Tint is assumed to be smooth and variable and the outer boundary r out Lipschitz 
and fixed. An example of such a domain is illustrated in figure [Tl] 

Remark 1 . In the shape optimization method proposed in this work the topology 
of the domain 12 is fixed. This means we do not consider topology optimization. 


Ttop 



Fig. 2.1: Example of a domain 12 = 12 1 U r^t U 12 2 where r out := <912 = Fbottom LI 
rleft U r right U rtop and n denotes the unit outer normal to S2 2 at T; nt 


The parabolic PDE constrained shape optimization problem is given in strong 
form by 


min J(f2) 

— j(^)+jreg(f2) : — 

-- (y — y) 2 dxdt + p Ids 

(2.6) 

dy 

sx - i - 



J 0 J Q J r i n t 


div(fcVy) = 

ymodel 

in 12 x (0, T\ 

(2.7) 


y = 

1 on 

r t o P x (0, T\ 

(2.8) 


dy _ 
dn 

0 on 

(Tbottom U Tieft U F r ig ht ) X (0, T\ 

(2.9) 


y = 

2/o in 

12 x {0} 

(2.10) 


where 


k = 


k\ = const. in 12f x (0, T] 
&2 = const. in 12 2 x (0, T] 


and n denotes the unit outer normal to 122 at Fj nt . Of course, the formulation (2.7) 


of the differential equation is to be understood only formally because of the jumping 
coefficient k. We observe that the unit outer normal to 12f is equal to —n, which 
enables us to use only one normal n for the subsequent discussions. Furthermore, we 
have interface conditions at the interface T; nt . We formulate explicitly the continuity 
of the state and of the flux at the boundary as 


= 0 , 


k ?y 

dn 


= 0 on r int x (0, T] 


( 2 . 11 ) 
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where the jump symbol [•] denotes the discontinuity across the interface Ti nt and 
is defined by [v] := V\ — where v\ := v\ and V 2 := v\ . The perimeter 

regularization, j leg (fl) := y f r t Ids, with p > 0 in the objective (2.6) is frequently 
used in this kind of problems. In |2&| a weaker but more complicated regularization 
is instrumental in order to show existence of solutions. 

We assume / model g L 2 (0,T; L 2 (fl)) and y € L 2 (0, T; Id 1 (fl)). In our setting, 
the boundary value problem (2.7|2.11) is written in weak form 


„ model / 


(y,p) = b™ de \p, P V) , v P g w (0 ,T-,H\n)) 


( 2 . 12 ) 


and for all p 1 € W (0, T; H 1 / 2 ^^ U T left U T right )), p 2 G W (0, T; Id- 1 / 2 (r top )) 


as m |^4j. For properties of the function spaces, we refer the r eader to the lit¬ 
erature, e.g. PU ESI- The bilinear form a mode> (y,p) in (|2.12[) is achieved by 


applying integration by parts on f Q ^j-p dx dt and on J Q f n div(kVy)p dx dt = 


So In, ■ 


cliv(fc 1 Vyi)pi dxdt + f Q f n div(fc 2 V?/ 2 )P 2 dxdt. Thus, we get 


model 


(y,p) ■= [ y(T,x)p(T,x)dx - [ y o p(0, x) dx - [ [ j^ydxdt 

Jn Jn Jo Jn at 


/0 Jn 
fT 


/0 J r D 


fcVr/ \7pdxdt — 

k\ ^-p dsdt. 
an 


,T , 

’ dy 

to J r int 

k^~P 

on 


ds dt 


(2.13) 


The linear form & model (p,p 1 ,p 2 ) in (2.12) is given by 


6 model (p,p 1 ,p 2 ) := b™° de \p) + br de \p\p 2 ) 


(2.14) 


where 



P° del pdxdt , 
p l (y — 1) ds dt + 





2 d y 


t\r top 


p —— ds dt. 
on 


(2.15) 

(2.16) 


In the following, we assume for the observation y 
grangian of (2.6||2.11) is defined as 


L 2 (0, T;iJ 1 (fl)). The La- 


-mj/,p) := JP) + a model (p,p) - b™ de \p,p\p 2 ) (2.17) 


where J(f2) is defined in (2.6), a model (y,p) in (2.13) and 6 model (p,p 1 ,p 2 ) in (2.14 


2.16). 


The adjoint problem, which we obtain from differentiating the Lagrangian 
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with respect to y, is given in strong form by 


dp 

dt 


div(fcVp) = 
V = 



dp 
dn 
V = 



{y-y) in 12 x [0, T) 

(2.18) 

in 12 x {T} 

(2.19) 

on Tint x [0, T) 

(2.20) 

on Tint x [0, T) 

(2.21) 

on (r bottom U Tjeft U r rig ht) x [0, T) 

(2.22) 

on T top x [0, T) 

(2.23) 

k\p on (T bo ttom C Tieft U T r i g ht) x [0, T) 

(2.24) 

dv 

L— on r top x [0, T) 

(2.25) 


and the state equation, which we get by differentiating the Lagrangian Jz? with 
respect to p , is given in strong form by 


dy 

dt 


div(fcVy) = / model in 12 x (0, T\. 


(2.26) 


As mentioned earlier, in many cases, the shape derivative arises in two equiva¬ 
lent forms. If we consider the objective (2.6) without the perimeter regularization 
jreg) the shape derivative can be expressed as an integral over the domain 1I as well 
as an inte gral over the interface Ti nt . Assume that a solution y of the parabolic PDE 
problem ( 2.7|2.11 ) exists and is at least in L 2 (0, T ; H 1 (12)). Moreover, assume that 
the adjoint equation (2. 18|2.23 ) admits a solution p G W (0,T; i7 1 (f2)). Then the 
shape derivative of the objective J without perimeter regularization, i.e., the shape 
derivative of j at 12 in the direction V expressed as an integral over the domain 12 
is given by 


Djn[V\ = f f - kVy T (W + VC T ) Vp - p ( V / model ) T 1/ 

Jo Jn 

+ div(V) (hy - y) 2 + ^p + fcVp T Vp - / model p ) dxdt. 


(2.27) 


This domain integral allows us to calculate the boundary expression of the shape 
derivative, which is given by 

Djr iat [V}= [ [ [ft] S7yJVp 2 ( V., n) ds dt. (2.28) 

Jo JT int 

The derivations are very technical. Note that we need a higher regularity of y and 


p to provide the boundary shape derivative expression (2.28). More precisely, p has 


to be an L 2 (0, T; if 2 ($2))-function having weak first derivatives in L 2 (0,T; H l { 12)') 
and y has to be an ele ment o f L 2 (0,T; H 2 (Q)). Here iJ 1 (f2)' denotes the dual space 
of i7 1 (f2). We achieve (2.27) by an application of the theorem of Correa and Seeger 
0 theorem 2.1] and (2.281 by an application of integration by parts. We refer the 
reader for its derivations to [24]. By combining theorem 2.1 and 2.2 in [24. with 
proposition 5.1 in m we get the following two expressions for the shape derivative 
of the objective J (with perimeter regularization) at 12 in the direction V: 


Djn [V) + Dj„ g (fi)[V] = DJ(Q)[V] = Dj Ti jV] + Dj n g (12)[H] 


(2.29) 
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with 


Dj res (Cl)[V] = / (V, n) g,n ds 


(2.30) 


/r in 


where n denotes the curvature corresponding to the normal n. 

3. Steklov-Poincare type metrics on shape manifolds. We first discuss 
the definition of shape manifolds and metrics. Then, we introduce novel metrics 
dovetailed to shape optimization based on domain formulations of shape derivatives. 

3.1. Shape manifolds. As pointed out in [23], shape optimization can be 
viewed as optimization on Riemannian shape manifolds and resulting optimization 
methods can be constructed and analyzed within this framework, which combines 
algorithmic ideas from pQ with the differential geometric point of view established 
in [16]. As in [23], we first study connected and compact subsets Q 2 C C R 2 with 
0 and C°° boundary where fl denotes a bounded domain with Lipschitz- 


boundary (cf. figure 2.11. We now identify the variable boundary = T; nt with a 
simple closed curve c: S l —> R 2 . Additionally, we need to describe a space including 
all feasible shapes Ti nt and the corresponding tangent spaces. In nsi, this set of 
smooth boundary curves c is characterized by 


B e (S\R 2 ) :=Emb(5 ,1 ,R 2 )/Diff(5' 1 ), 


(3.1) 


i.e., as the set of all equivalence classes of C°° embeddings of S 1 into the plane 
(Emb)^ 1 , R 2 )), where the equivalence relation is defined by the set of all C°° re- 
parameterizations, i.e., diffeomorphisms of S 1 into itself (Dif^S 1 )). A particular 
point on the manifold B e (S 1 ,'R 2 ) is represented by a curve c: S 1 9 6 >->■ c(6 ) £ R 2 . 
Because of the equivalence relation (Diff(5 1 )), the tangent space is isomorphic to 
the set of all normal C°° vector fields along c, i.e., 

T c J5 e = {h: h = an, a £ C°°(S' 1 ,R)} (3.2) 


where n is the unit exterior normal field of the shape defined by the boundary 
dfl 2 = c such that n(9) T c' for all 6 £ S 1 and d denotes the circumferential 
derivative as in m- Several intrinsic metrics are discussed in m , among which 
the following Sobolev metric seems the most natural intrinsic one from a numerical 
point of view. For A > 0, the Sobolev metric is induced by the scalar product 

g 1 '■ T c B e x T c B e —> R, 

( h , k) h> ((id - AA c )a, p) L 2 (c) 

where h = an and k = /3n denote two elements from the tangent space at c and A c 
denotes the Laplace-Beltrami operator on the surface c. In [16( it is shown that the 
condition A > 0 guarantees that the scalar product g 1 defines a Riemannian metric 
on B e and thus, geodesics can be used to measure distances. 

With the shape space B e and its tangent space in hand we can now form the 
Riemannian shape gradient corresponding to a shape derivative given in the form 


dJ(fi)[E] = 


ip (V, n) ds. 


(3.4) 


In our model setting the objective function J is given in (2.6) and its shape deriva¬ 
tive in (2.29). Finally, the Riemannian shape gradient gradJ with respect to the 
Riemannian metric g 1 is obtained by 


gradJ = gn with (id — AA c )g = ip . 


(3.5) 
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The metric g 1 , which is also used in [24] . necessitates a shape derivative in Hadamard 
form as an efficient means to solve linear systems involving the Laplace Beltrami 
operator in surfaces. All of that is certainly not impossible but requires computa¬ 
tional overhead which we can get rid of by usage of the metric discussed below. We 
compare the algorithmic aspects of both approaches below in section [5] 

3.2. Steklov-Poincare type Riemannian metrics. The ideal Riemannian 
metric for shape manifolds in the context of PDE constrained shape optimization 
problems is to be derived from a symmetric representation of the second shape 
derivative in the solution of the optimization problems. Often, this operator can 
be related to the Dirichlet to Neumann map, aka Steklov-Poincare operator, or 
the Laplace-Beltrami operator [21] , If one aims at mesh independent convergence 
properties, one of these two will be appropriate in most cases. Since it can be 
observed that the Laplace-Beltrami operator is spectrally equivalent to the square 
of the Steklov-Poincare operator, the latter operator seems to be more fundamental 
and we will focus on it as a basis for the scalar product on T c B e . Another advantage 
of this operator is that is blends well in with a corresponding mesh deformation 
strategy. 

Most often, the Dirichlet to Neumann map is associated with the Laplace op¬ 
erator. However, as pointed out in mm more general elliptic operators can be 
involved. For the purpose of mesh deformation, an elasticity operator may be the 
ideal choice. In numerical computations, its inverse, the Neumann to Dirichlet map 
or Poincare-Steklov is also of importance. Therefore, we first define these operators. 

In the sequel, we use the continuous generalized trace map 


7 : H^{Sl,R d ) H 1 / 2 (T int ,M d ) x R-V 2 (r int ,R d ), 


U 1 —^ 


f 7oEA 

\rnu) 



(3.6) 


Analogously to IL3], we define for vector fields I/, V £ (Q, M d ) with d = 2 or 

d = 3, the Neumann solution operator for the inner boundary Ti nt derived from a 
symmetric and coercive bilinear form 


by 


a: H^(Q,R d ) x H^(Q,R d ) -> ffi. (3.7) 

E n : iL _ 1 // 2 (r int ,R d ) -> H^{Sl,R d ), 


with U defined as the solution of the variational problem 


a(U,V) = [ u T { 70 F) ds, V V G Hg(Q,R d ) 


(3.9) 


where we note that the integral in the right hand side of equation (3.91 is to be 
understood as the duality pairing. Furthermore, we define the Dirichlet solution 
operator for the inner boundary r; nt by 


E d -. H 1 ^(T illt ,R d )-^H^n,R d ), 
u 1 —^ U 

with U defined as the solution of the variational problem 

a (U,V) = 0, Vkelo^,!'); u\ 

' -L int 


(3.10) 


= u. 


(3.11) 



Now, we can define the Dirichlet to Neumann map and the Neumann to Dirich- 
let map as done in the following definition: 

Definition 3.1. In the setting above, the Dirichlet to Neumann map T and 
the Neumann to Dirichlet map S are defined by 


T := 7l o E d : H 1 ^ 2 (r int ,'R d ) ITV 2 (r int ,I 
5 := 70 o E n : ff- 1 / 2 (r int , K d ) i? 1 / 2 (r int ,I 


*), 


(3.12) 

(3.13) 


where 7o,7i are given in (3.6). 


In obvious generalization of theorem 2.3.1 in |l3j from scalar fields to vector 
fields, we conclude that both operators are symmetric w.r.t. the standard dual 
pairing (■,■), coercive, continuous and that T = S~ x , an observation, for which [2] 
is cited in m ■ For the purpose of defining an appropriate scalar product on the 
tangent space of shape manifolds, we define the following mappings. 

Definition 3.2. In the setting above, we define 


ry. H(T iat ) —► H(T inU 
a ia ■ n 


*) V T - H(T int ,R d ) ^ H(T iDt ) 


U i-> n T U 


where H £ {H 1//2 , if 1 / 2 }, and thus the projected operators 

TP := V T o T o ry. ff 1 / 2 (r int ) -> I?- 1 / 2 (r int ), 
S p := r, T o S o ry. H~ l ' 2 (T- mt ) H l / 2 (V int ). 


(3.14) 

(3.15) 


Both operators, T p and S p , inherit symmetry, coercivity, continuity and invert- 
ibility from the operators T, S. However, we observe in general T p ^ ( S p ) _1 . Both 
operators can be used for the definition of a scalar product on the tangent space. 
In line with the discussion of Sobolev type metrics in [16] , we would prefer a scalar 
product with a smoothing effect like the projected Dirichlet to Neumann map T p . 
However, we need its inverse in numerical computations, which is usually not S p , 
although spectrally equivalent. We can limit the computational burden, if we use 
directly (5 P ) _1 as a metric on the tangent space, having a similar smoothing effect 
but also the advantage of the straight forward inverse S p . In order to summarize, 
let us explicitly formulate the operator 


(3.16) 


S p : # _ 1 / 2 (r int ) ->• H^ 2 (Ti nt ), 

a i ^ (yoU) T n 
where U € Ffo(D,M d ) solves the Neumann problem 

a(U, V) = f a-^ 0 V) T n ds, V V £ H({n,R d ) (3.17) 

which corresponds to an elliptic problem with fixed outer boundary and forces a ■ n 
at the inner boundary Tint. Thus, we propose to use the scalar product g s defined 
below. 

Definition 3.3. In the setting above, we define the scalar product g s on 
ff 1/2 (r int ) by 

g s : iJ 1 / 2 (ri nt ) x iJ 1 / 2 (r int ) — ► R, 

(a,p) < a ,(S*’)-V) = [ a(s) ■ [(S p )-^](s) ds. 

•I r int 


(3.18) 
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4. Shape quasi-Newton methods based on the metric g s . As already 
mentioned in section [2] the shape derivative can always be expressed as boundary 
integral DJr iat [V] = / r . ^f{V,n)ds (cf. (2.4)) due to the Hadamard structure 


theorem. If V | ^ = an, this can be written more concisely as 


DJr iat [V}= f afds. 


(4.1) 


Due to the isomorphism ( |3.2| ) and the handy expression (4.1) we can state the 
connection of (BgfS 1 , R 2 ), g a ) with shape calculus, i.e., we can determine a repre¬ 
sentation h £ Tr^BejS 1 , R 2 ) = {h: h = an, a £ C°°(S 1 , R)} of the shape gradient 


(4.2) 


in terms of g s defined in (3.18) by 

9 S {^h) = (/,</>) L 2 (rint) , V</> £ C°° (r int ,] 

which is equivalent to 


(4.3) 


[ fi(s) ■ [(S ,p )- 1 /i](s) ds= [ f 000(a) da, V0 £ C°°(r int ,R). 
dr int J r int 

Thus, h = S p f = ( 7 oD) T n, where U £ 1^(12, R d ) solves 
a(U, V) = f / • ( 7o V) T n ds = DJ Fi jV\ = DJ n [V} ,Vb£ (4.4) 

J rint. 


which means that the representation of the domain integral formulation in terms of 
the elliptic form a(-, •) as used in m can be - projected to the normal component 
on Tint - interpreted as the representation of the boundary integral formulation in 
terms of (5 lp ) _1 . However, in both points of view, the information of the shape 
derivative is in physical terms used as a force (in the domain or on the boundary) 
and we obtain a vector field U as an (intermediate) result, which can serve as a 
deformation of the computational mesh identical to Dirichlet deformation. 

Remark 2. In general, h = S p f = (7o U) T n is not necessarily an element of 
Tr int B e because it is not ensured that U £ .ffg(n,M d ) is C°°. Under special assump¬ 
tions depending on the coefficients of a second-order partial differential operator, the 
right hand-side of a PDE and the domain fl on which the PDE is defined, a weak 
solution U £ .ffg(fl,R d ) of a PDE is C°° by the theorem of infinite differentiability 
up to the boundary theorem 6, section 6.3]. 

Now, we rephrase the 1-BFGS-quasi-Newton method for shape optimization 
from (241 in terms of the metric g s and in generalization to domain formulations 
of the shape derivative. We note that the complete deformation of a shape opti¬ 
mization algorithms is just the (linear) sum of all iterations, which means that the 
BFGS update formulas can be rephrased directly in terms of the deformation vector 
field, rather than only as boundary deformations to be transferred to the domain 
mesh in each iteration. 

BFGS update formulas need the evaluation of scalar products, where at least 
one argument is a gradient-type vector. According to the metric introduced in 
section |3j we can assume that a gradient type vector u £ T c B e can be written as 

u = (7o U) T n (4.5) 


for some vector field U £ The other argument v is either of gradient- 

type or deformation-type, which can also be assumed of being of the form (4.51, 
i.e., 


v = (jo V) T n 


(4.6) 
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for some V £ i?o(f2,M d ). If u is a gradient of a shape objective J, we observe 

/(«,«) = DJ Tini [V] = DJ n [V ] = a(t/,R). (4.7) 

This observation can be used to reformulate the scalar product g s (■,■) on the bound¬ 
ary equivalently as a(-, •) for domain representations. In the sequel, we consider only 
domain representations Uj £ I7g(fi,IR d ) of gradJ(cj) £ H 1 / 2 { r int ), mesh deforma¬ 
tions Sj £ 77Q(fl,K d ) and differences Yj := Uj + 1 — Ts^Uj £ .Hg(f2,M d ) where Ts t 
denotes the vector transport as in [24]. 

With this notation we formulate the double-loop of an 1-BFGS quasi-Newton 
method: 

Pj 9 S (hoYj) T n, ( 70 Sj) T n) _1 = a(Y), S,-) -1 
q £- I/,- 

for i = j — 1 ,,j — m do 

Si <- r g s, 

<- T q Yi 

a* £- /^g 5 (( 7 oS’i) T n, (7og) T R-) = p*o(Si, g) 
q «— q - onYi 

end for 

g S (hoYj-i) T n,( r yoS j -i) T n) _ a(yj-i,Sj-i) 

^ S s ((7o>'j-i) T ™,(7oy,--i) T >T.) ^ _ a(Yj -” 

for i = j — m,... ,j — 1 do 

Pi <- P*5 S ((7o^) T? b ( 7 o^) T n) = pia(Yi,q) 
q <- q + {at - Pi)Yi 

end for 

return g = Gj 1 gradJ(c J ) 

The resulting vector g £ Ffg(f2,R d ) is simultaneously a shape deformation as well 
as a deformation of the domain mesh. 


5. Numerical results and implementation details. We compare the lim¬ 
ited memory BFGS shape optimization algorithms of [24] with the analogous algo¬ 
rithm based on the Riemannian metric g s , introduced above. We use a test case 
within the domain f l = (—1, l) 2 , which contains a compact and closed subset 1^2 
with smooth boundary. The parameter k\ is valid in the exterior fR = Q \ Q 2 
and the parameter k 2 is valid in the interior 0,2- First, we build artificial data y, 
by solving the state equation for the setting := { x : ||x|| 2 < r} with r = 0.5. 


Afterwards, we choose another initial domain IR and fl 2 - Figure [2T] illustrates the 
interior boundary around the initial domain fl 2 and the target domain The 
reason for this choice of artificial test data is that we obtain a representation of y 
that can be evaluated at arbitrary points in space since it is represented in finite 
element basis. Moreover, this construction guarantees that the optimization con¬ 
verges to a reasonable shape that is within the boundaries Q and not too different 
to the initial shape such that the mesh remains feasible under deformations. 

Remark 3. Choosing the measurements y as the solution of the model equation 
{'2. 1 2.1l\) we obtain that y £ L 2 (0,T; 17 1 (f2)) as we assumed in section 


2.1 


In the particular test case, which is studied in this section and can be seen in 
figure [573[ the diffusion coefficients are chosen to be £q = 1 and ^2 = 0.001. Further, 
the initial condition is yo(x) = 0 for all x £ Cl, f nwdel (x,t ) = 0 in (x, t) £ fl x (0, T] 
and the final time of the simulation is T = 20. The results shown in this section 
are computed under a mild perimeter regularization with y = 10 -6 , where we did 
not notice any numerical difference with the case y = 0. Yet, for the non-smooth 
initial configuration shown in figure |5.4| a stronger regularization has to be chosen 
in the first iterations as /q n it = 0.01. In this particular case the regularization is 
controlled by a decreasing sequence from /q n i t to y. 

The numerical solution of the boundary value problem ( 2.7||2.10 ) is obtained by 
discretizing its weak formulation (2.12) with linear finite elements in space and an 
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implicit Euler scheme in time. For the time discretization 30 time steps are chosen, 
which are equidistantly distributed. The diffusion parameter k is discretized as a 
piecewise constant function in contrast to the continuous trial and test functions. 


This choice of function spaces ensures that the transmission conditions (2.111 are 
automatically fulfilled. The corresponding adjoint problem ( 2.18||2.25 ) can be dis¬ 
cretized in the same way. More precisely, it is not necessary to assemble different 
linear operators, which is attractive in terms of computational effort. All arising 
linear systems are then solved using the preconditioned conjugate gradient method. 

Our investigations focus on the comparison between two 1-BFGS optimization 
approaches: The first approach is based on the surface expression of the shape 
derivative, as intensively described in [24]. Here, a representation of the shape 


gradient at Fj nt with respect to the Sobolev metric (3.3) is computed and applied 


as a Dirichlet boundary condition in the linear elasticity mesh deformation. This 
involves two operations, which are non-standard in finite element tools and thus 
leads to additional coding effort. Since we are dealing with linear finite elements 
the gradient expressions of state y and adjoint p in ( |2.28 l are piecewise constant 
and can not be applied directly to the mesh as deformations. We thus have to 
implement a kind of L 2 -projection on r; nt (cf. [M]) bringing back the sensitivity 
information into the space of continuous, linear functions. The next additional piece 
of code is a discrete version of the Laplace-Beltrami operator for the Sobolev metric 


(3.3). The essential part of this is the solution of a tangential Laplace equation on 
the surface Tint. Therefore, we follow the presentations [T5] and artificially extend 
our 2D grid in the third coordinate direction. The second approach, discussed 
in sections [3] and [4j involves the volume formulation of the shape derivative and a 
corresponding metric, which is very attractive from a computational point of view. 
The computation of a representation of the shape gradient with respect to the chosen 
inner product of the tangent space is now moved into the mesh deformation itself. 
The elliptic operator a(-,-) (cf. (3.7)) here the linear elasticity - is both used as 
inner product and mesh deformation leading to only one linear system, which has to 
be solved. Besides saving brain work in the calculation of the shape derivative, a lot 
of coding work is obsolete using surface formulation of shape derivatives. Moreover, 
it is not always clear how the surface formulation looks like and which additional 
assumptions have to be made in its derivation. A discussion of the 1-BFGS algorithm 
used within this algorithm can be seen in section [4] 

An essential part of a shape optimization algorithm is to update the finite 
element mesh after each iteration. For this purpose, we use a solution of the linear 
elasticity equation 


div(cr) = / olas in 
U = 0 on F out 

a := ATr(e)/ + 2pe 

e:=l(VU + VU T ) 


(5.1) 


where a and e are the strain and stress tensor, respectively. Here A and p denote 
the Lame parameters, which can be expressed in terms of Young’s modulus E and 
Poisson’s ratio v as 


A = 


vE 


(l + i/)(l — 2 v) 


M = 


E 

2(1 + v) 


(5.2) 


The solution U \ SI —¥ R 3 is then added to the coordinates of the finite element 
nodes. The Lame parameters do not need to have a physical meaning here. It is 
rather essential to understand their effect on the mesh deformation. E states the 
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stiffness of the material, which enables to control the step size for the shape update. 
v gives the ratio how much the mesh expands in the remaining coordinate directions 
when compressed in one particular direction. The numerical results in this work 
are obtained using v — 0.01 and E = 0.1. 


Equation (5.1) is modified according to the optimization approach under con¬ 
sideration. In case we use the surface formulation of the shape derivative (2.28), 
the following Dirichlet condition is added on the variable boundary 


U = U 


surf 


on 


r ir 


(5.3) 


where U suli is the representation of the shape gradient with respect to the Sobolev 


metric g 1 given in (3.3). The source term / Glas is then set to zero. Otherwise 


when the mesh deformation operator is also used as shape metric, / elas assembled 


according to (2.27) and there is no Dirichlet condition on U. This only covers the 


portion of the shape derivative for which a volume formulation is available. Parts of 
the objective function leading only to surface expressions, such as, for instance, the 
perimeter regularization j reg , are incorporated as Neumann boundary conditions 
given by 


^^ rsurf 

dn ~ 1 


on 


int* 


(5.4) 


In the notation of section 3.2 we set ad¬ 
equation leading to 


•) as the weak form of the linear elasticity 


a(U, V) = f a{U) : e(V) dx. (5.5) 

J n 

For our model problem given in section [2] we have to solve in the context of the 
domain formulation of the shape derivative and its representation in terms of g s 

a(U , V) = Dj n [V) + Dj ieg (n)[V] , V V G H*(n,R d ) (5.6) 


where the right hand side is given by the left formulation in (2.29), in particular 
Djn[V] in (2.27) and Dj Ieg (Cl)[V] in (2.30). Note that (5.6) is justified by the main 
result (4.4) of the previous sections stating that the connection between the volume 
formulation of shape derivatives and a bilinear form a leads to a representation of 
the shape gradient with respect to g s . 

Both approaches (A versus B below) follow roughly the same steps with a major 
difference in the way the shape sensitivity is incorporated into the mesh deformation. 
The appraoch A (domain formulation) is clearly to be preferred because of its 
implementational ease and less computational effort, if a technical detail discussed 
below is taken into account. One optimization iteration can be summarized as 
follows: 


1. The measured data y(t,x) has to be interpolated to the current iterated 
mesh and the corresponding finite element space. Here, this consists of the 
interpolation between two finite element spaces on non-matching grids. 

2. The state and the adjoint equation are solved. 

3. Assembly of the linear elasticity equation. 

A) Domain formulation: 

• The volume form of the shape derivative is assembled into a source 
term for the linear elasticity mesh deformation. Only test functions 
whose support includes Ti nt are considered, which is justified in the 
subsequent discussion. The behavior of the algorithm with full as¬ 
sembly for all test functions is illustrated in figure |5.2b[ Here, the 
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magnitude of the unmodified discretization of the source term is vi¬ 
sualized, which shows not only non-zero values outside of Fi nt due to 
discretization errors, but leads also to detrimental mesh deformations. 

• Shape derivative contributions, which are only available in surface 
formulation, such as the perimeter regularization, are assembled into 
the right hand side in form of Neumann boundary conditions. 


B) Surface formulation: 


• The preliminary gradient s = Vy-[ Vpyn given in (2.28) is evaluated 

Tj n t. 

• The L 2 -projection of s into the space of piecewise linear, continuous 
functions is conducted. Let this be denoted by s. 

• Finally the contributions resulting from the regularization, which is 
here «n, is added to s and we solve the Laplace-Beltrami equation 
(id — AA c )U surt = s + nn to obtain a representation of the gradient 
with respect to the Sobolev metric as given in (3.3). 

• [7 surf then yields the Dirichlet boundary condition (5.3). 

4. Solve linear elasticity equations, apply the resulting deformation to the 
current finite element mesh and go to the next iteration. 


Assembling the right hand side of the discretized weak form (equation ( |5.6| )) 
only for test functions whose support intersects with T; nt in the volume formulation 
of step 3 above is due to the following reasoning. In exact integration, the integral 
Djn [V] should be zero for all test functions V which do not have Ti nt within their 
support. Thus, nonzero integral contributions are caused by discretization noise, 
On the other hand, its effect on the optimization algorithm can be understood from 
a perturbation point of view. We may assume that the Riemannian shape Hes¬ 
sian V c gradJ (where V c means covariant derivative), whose action in the optimal 
solution coincides with the action of the shape Hessian, i.e., 


/(V c gradJ[R], U) = D(DJ[V})[U] (5.7) 


is coercive on the boundary, i.e., for projections 77 T R|r int , ^ T H|r int , which guar¬ 
antees a well-posed problem. However, the Hessian operator approximated in the 
BFGS update strategy described in section [4] deals with a Hessian defined on the 
whole mesh, which posseses a huge kernel, determined by all vector fields with zero 
normal component on the boundary. Thus, the space ffQ(fi,]R d ) of all admissible 
deformations has a decomposition 


H'(n,R d ) = 


Hy 

L in 


(5.8) 


where Hy iiit := {£W(cm): a € iJ -1 / 2 (Tint)} and denotes its orthogonal com¬ 

plement in the bilinear form a(-,-). Shape gradients and increments in iJQ(H,K d ) 
lie in Hr int only. It is abvious that 1-BFGS update formulas produce steps which lie 
again in Hy. nl only, which means that the optimization algorithm in function spaces 
acts always on the coercive shape Hessian only. However, the discretized version 
is a perturbation of the infinite Hessian. Thus, perturbed coercive operators stay 
coercive, if the perturbation is not too large. But, positive semidefinite operators 
with a nontrivial kernel, almost inevitably will get directions of negative curvature, 
when perturbed. These directions of negative curvature will be chosen, if we al¬ 
low nonzero components in the right hand side of the discretized mesh deformation 


equation (5.6) in the interior of the domain. On the other hand, if we do not allow 


zero components there, the algorithm only acts in the subspace of the discretization 
of Hy iat where the projected Hessian is a perturbation of the shape Hessian and 
thus coercive, if the perturbation is not too large. 
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Fig. 5.1: Convergence of the optimization iteration measured as an approximation 
to the geodesic distance in the shape space on a grid with approx. 100,000 cells 



(a) BFGS iterates with unmodified approximation of volume shape derivative 



(b) Magnitude of unmodified volumic source term 


Fig. 5.2: Wrong mesh deformations and source term due to discretization errors in 
the unmodified right hand side of (5.6 1 


We now conclude this section with a discussion of the numerical results. The 
figures |5.2| to |5.4| show the initial configuration and the iterations 2, 4, 20 of the 
full BFGS algorithm as described in section |4j In figure |5.2| the algorithm is shown 
for the unmodified assembly of the right hand side in ( |5.6[ ) leading to divergence. 
Whereas, figure |5.3| shows a selection of BFGS iterates for the modified source 
term and with fj, = 0. It is also demonstrated here that the domain-based shape 
optimization algorithm can be applied to very coarse meshes. This is due to the 
fact, that there is no dependence on normal vectors like in the case of surface shape 
derivatives. Finally, figure p)~T| shows the convergence of 1-BFGS with three gradients 
in memory, full BFGS and the pure gradient method for suface and volume shape 
derivative formulation, respectively. 
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• 



’a) Smooth deformations and convergence to 

optimal shape due to modified source term 

* 


• 


• 


• 


(b) The domain-based optimization approach also enables the use of much coarser spatial 
discretizations (approx. 1000 cells) 


Fig. 5.3: BFGS iterates with corrected source term (2.27) indicating mesh indepen¬ 
dent convergence 






Fig. 5.4: Smooth mesh deformations even with kinks in the initial configuration due 
to regularization 


In our tests, the convergence with the Laplace-Beltrami representation of the 
shape gradient seems to require a bit less iterations compared to the domain-based 
formulation. Yet, the domain-based form is computationally more attractive since 
it also works for much coarser discretizations. This can be seen in figure [Y3| where 
|5.3a| shows the necessary fineness of the mesh for the surface derivative to lead to 
a reasonable convergence. The coarse grid in 5.3b however, only works for the 
domain-based formulation. 

Since the volume term -DjnjV] +Dj leg (fl)[V] in approach A only has to be com¬ 
puted for discretization elements adjacent to Ti nt it is computationally not more 
expensive than the surface formulation in approach B. Moreover, the computing 
time for L 2 -projection and solution of tangential Laplace equation is saved. Yet, 
the BFGS update algorithm in approach A is more expensive then the one in B due 
to the higher dimension of the involved matrix, which does not play a decisive role 
since we only store a few gradients in memory. These differences should yet not be 
overrated. The most expensive operation is the computation of the mesh deforma¬ 
tion involving the solution of the linear elasticity equation in both approaches A 
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and B making them comparable in terms of computational costs. 

This changes for highly parallel application on supercomputers as investigated 
in HZ]. Operations, which are only performed on surfaces, can drastically affect the 
scalability of the overall algorithm if the computational load is not balanced also 
with respect to surface elements. The higher demand for memory of the domain- 
based formulation seems also not be dramatic since the numerical tests suggest that 
very few gradients in memory are sufficient for good performance of the 1-BFGS 
method (see figure |5T| ). 

6. Towards a novel shape space. The scalar product introduced above 
in section [ 3 ] connects shape gradients with H 1 deformations. These deformations 
evaluated at a prior shape T 0 give deformed shapes Tj nt of class H 1 / 2 , if the defor¬ 
mations are injective and continuous. In the following, it is clarified what we mean 
by 77 1 / 2 -shapes. The investigations done in the previous section are not limited to 
C°° shapes, i.e., elements of the shape space i3 e (5' 1 ,R 2 ). Therefore, this section is 
devoted to an extension of B e (S 1 , R 2 ), i.e., to a novel shape space definition, and 
its connection to shape calculus. 

We would like to recall once again that a shape in the sense of the shape space 
of Peter W. Michor and David Mumford introduced in [16] is given by the image of 
an embedding from the unit sphere S^ -1 into the Euclidean space R d . In view of 
our generalization, it has technical advantages to consider a prior shape T 0 as the 
boundary T 0 = dX 0 of a connected and compact subset X 0 C O C R d with X 0 ^ 0, 
where Q denotes a bounded Lipschitz domain. Let the prior set X$ be a Lipsclritz 
domain, i.e., To is a Lipschitz boundary. An example of a prior shape is the cube. 
It is the union of six faces, where each is a portion of a plane, i.e., a smooth surface. 
General shapes - in our novel terminology - arise from if 1 -deformations of such a 
prior set X 0 . These iJ 1 -deformations, evaluated at a prior shape T 0 = dX 0 , give 
deformed shapes Ti nt if the deformations are injective and continuous. We call these 
shapes of class H 1 / 2 and define the set 

-H 1/2 (T 0 ,R d ) := {w: T 0 -» A': 3W <E H\n,n) s.t. 

W\ injective, continuous, W\ = w}. 
r 0 T 0 


However, in order to have a unique representation for each shape, we have to factor 
out the homeomorphisms from the prior shape T 0 into itself which are compatible 
with the set (6.1). Thus, we characterize the following shape space: 

Definition 6.1. Let Q, Xq and To be as above. The space of all H 1 ^ 2 -shapes 
is given by 


H 1/2 (T 0 ,R d ) := H 1 / 2 (r 0 ,R d )/Homeo 1/2 (r 0 ), 


( 6 . 2 ) 


where 'H 1 / 2 (ro,R d ) is given in \(>. l\ ) and Homeo 1 ^ 2 (To) is defined by 

Homeo 1 ^ 2 (T 0 ) := {ru: w £ 'H 1 ^ 2 (T 0 ,R d ), w: T 0 —» T 0 homeomorphism}. (6.3) 


Remark 4. Of course, the properties of the shape space B 1 ^ 2 (To,R d ) have to be 
investigated. For example the independence of the prior shape T 0 in the shape space 
definition is an open question. If it is independent, we can choose, for example, the 
unit sphere S d ~ 1 as prior shape. Another important question is whether the shape 
space has a manifold structure. Note that this question is very hard and a lot of 
effort has to be put into it to find the answer. From a theoretical point of view there 
are several other open questions. However, this goes beyond the scope of this work 
and is a topic of subsequent work. 

Remark 5. In the following, we assume that B 1 ^ 2 (To,R d ) has a manifold 
structure. If necessary, we can refine the space S 1 / 2 (To,R d ), e.g., by restriction 
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to an explicit deformation field W. In our setting, it arises from the linear elas¬ 
ticity equation and the request of the existence of an arbitrary one is perhaps too 
strong. This way, we can replace by a linear space, which is in par¬ 

ticular a manifold. However, this conceivable limitation leaves the following theory 
untouched. 

if r e s 1 / 2 ^,!^) is smooth enough to admit a normal vector field n, the 


following isomorphisms naturally arise out of definition (6.1): 


T r B 1/2 (r 0 ,K d ) 

= {h: h = 4>n a.e., (j> G H 1 ^ 2 ( r) injective, continuous} 
= {<j>: <j> G H 1/2 (T) injective, continuous} 


(6.4) 


Now, we can formulate the shape quasi-Newton method s of section [4| on the 
novel shape space B 1 ^ 2 (r 0 ,R d ). Due to the isomorphism (6.4) and the handy 


expression (4.1) of the shape derivative we can state the connection of B 1 ^ 2 (F 0 ,1 
with respect to g s to shape calculus, i.e., we can determine a representation h G 
Tr int S 1 / 2 (r 0 ,R d ) = {h: h G H 1 ^ 2 ( r) injective, continuous} of the shape gradient 
in terms of g s defined in (3.18) by 


g s ((/),h) = (r,<j>) L 2 (r int) 


(6.5) 


for all injective and continuous <fi G H 1 / 2 ( Fi nt ), which is equivalent to 



</>(s) ■ [(S' p ) 1 h](s) ds 



r(s)<p(s) ds 


( 6 . 6 ) 


for all injective and continuous f> G H 1 ^ 2 ( ri nt ). 

Based on the connection ( |6.5[ ) we can f ormu late the quasi-Newton methods of 
section}} also on ( B x ! 2 (To,R d ) , g s ). From (6.6) we get h = S p r = ( 7 oU) T n where 
U G Hfflh,R d ) solves 


a(U, V)= r • ( 7 o V) T n ds = DJ Fiat [V] = DJ n [V ], VVG HZ( fi, R d ). (6.7) 
Jr int 

In general, h = S p f = ('yoU) T n is not necessarily an element of Tret'S 1 / 2 (r 0 ,R d ) 
because it is not ensured that U G i7Q(D,R d ) is injective and continuous. Under 
special assumptions depending on the coefficients of a second-order partial differ¬ 
ential operator, the right hand-side of a PDE, the domain f l on which the PDE is 
defined and the dimension of U, the continuity of a weak solution U G HQ(U,M. d ) 
of a PDE is guaranteed by the theorem of higher interior regularity [S] theorem 5, 
section 6.3] combined with the Sobolev embedding theorem. In particular, if these 
conditions are fulfilled, we get in our two- and three-dimensional case a bounded 
C 2 regularity of U. Now, if we require 

nic}(cW) < 1 (6-8) 

we get the injectivity of U due to 01 lemma 6.13 and remark 6.14]. 

Remark 6. In implementations, the necessary condition || CF||cj(o,]R d ) < 1 f or 
injectivity of the deformation U is ensured by the particular choice of the Lame 
parameters. 

7. Conclusions. This paper develops an intrinsic metric in shape spaces, 
which enables to jointly work with domain based and boundary based shape deriva¬ 
tive expressions, and which leads to shape optimization algorithms with several 
computational and analytic advantages as outlined above. Furthermore, the metric 
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leads to a novel shape space B 1 ^ 2 . The properties of B 1 ^ 2 are beyond the scope of 
this work and will be touched in subsequent papers. It is obvious that the results 
of this paper are not restricted to two space dimensions and also not to interior 
interface shapes. The whole discussion carries over to shapes which are just parts 
of the exterior boundary of a computational domain. 
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